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Abstract. - The calculation of the heating rate of cold atoms in vibrating traps requires a theory 
that goes beyond the Kubo linear response formulation. If a strong "quantum chaos" assumption 
does not hold, the analysis of transitions shows similarities with a percolation problem in energy 
space. We show how the texture and the sparsity of the perturbation matrix, as determined 
by the geometry of the system, dictate the result. An improved sparse random matrix model is 
introduced: it captures the essential ingredients of the problem, and leads to a generalized variable 
range hopping picture. 
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The rate of energy absorption by particles that are con- 
fined by vibrating walls was of interest in past studies of 
nuclear friction [1-3] , where it leads to the damping of the 
wall motion. More recently it has become of interest in 
the context of cold atoms physics. In a series of experi- 
ments [4-6] with "atom-optics billiards" some surprising 
predictions [7] based on linear response theory (LRT) have 
been verified. 

In the present study we consider the case where the bil- 
liard is fully chaotic [a] but with nearly integrable shape 
(Fig.l). We explain that in such circumstances LRT is 
not applicable (unless the driving is extremely weak such 
that relaxation dominates). Rather, the analysis that is 
relevant to the typical experimental conditions should go 
beyond LRT, and involve a "resistor network" picture of 
transitions in energy space, somewhat similar to a perco- 
lation problem. Consequently we predict that the rate of 
energy absorption would be suppressed by orders of mag- 
nitude, and provide some analytical estimates which are 
supported by a numerical calculation. 

We assume that an experimentalist has control over the 
position (R) of a wall element that confines the motion 
of cold atoms in an optical trap. We consider below the 
effect of low frequency noisy (non periodic) driving. This 
means that R is not strictly constant in time, either be- 
cause of drifts [8] that cannot be eliminated in realistic 
circumstances, or else deliberately as a way to probe the 
dynamics of the atoms inside the trap [9]. We assume 
the usual Markovian picture of FGR transitions between 



energy levels, which is applicable in typical circumstances 
(see e.g. [10]). These transitions lead to diffusion in the 
energy space. If the atomic cloud is characterized by a 
temperature T, then the diffusion in energy would lead 
to heating with the rate E = D/T [b] and hence to an 
increase in the temperature of the cloud. 

Naively one expects to observe an LRT behavior. That 
means to have D cx [RMS(i?)] 2 , and more specifically to 
have a linear relation between the diffusion coefficient and 
the power spectrum of the driving: 



D 



G x RMS(7?) 5 



C(u)S{uj)duj (1) 



Here S(uj) is the power spectrum of R, and C(u) is related 
to the susceptibility of the system. From the experimen- 
talist's point of view the second equality in Eq.(l) can be 
regarded as providing a practical definition for C(uj), if 
the response is indeed linear. 

We shall explain in this paper that the applicability 
of LRT in our problem is very limited, namely LRT 
would lead to wrong predictions in typical experimental 
circumstances. Rather we are going to use a more re- 
fined theory, which we call semi-linear response theory 
(SLRT) [11, 12], in order to determine D. The theory is 
called SLRT because on the one hand the power spectrum 
S(lli) h- > \S(ui) leads to D i— > XD, but on the other hand 
S(tu) i— > Si(uj) + S>2(u) does not lead to D ^ D\ + D 2 . 
This semi-linearity can be tested in an experiment in or- 
der to distinguish it from linear response. Accordingly, in 
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SLRT the spectral function C(u) of Eq.(l) becomes ill de- 
fined, while the coefficient G is still physically meaningful, 
and can be measured in an actual experiment. 

If we assume small driving amplitude the Hamiltonian 
matrix can be written as TL = {E n } + f{t){V nm }, where 



dU 



dR 



(2) 



is the perturbation matrix. More than 50 years ago 
Wigncr had proposed to regard the perturbation matrix 
of a complex system as a random matrix (RMT) whose 
elements are taken from a Gaussian distribution. Later 
Bohigas had conjectured that the same philosophy applies 
to quantized chaotic systems. For such matrices the va- 
lidity of LRT can be established on the basis of the FGR 
picture, and the expression for G is the Kubo formula 
G LRT = 7rg E ((|K l m| 2 ))a where {{x)) a = (x) is the algebraic 
average over the near diagonal matrix elements [c] , and g E 
is the density of states (DOS). In contrast to that, using 
the Pauli master equation [10] with FGR transition rates 
between levels, the SLRT analysis leads to 



= KQ E ((\V nm \ 2 )) 



(3) 



where the "average" ((x)) is defined as in Rcf. [11, 12] 
via a resistor-network calculation [13]. (For mathematical 
details see "the SLRT calculation" paragraph below) . 

Within the RMT framework an element x of | V^ rrl | 2 
is regarded as a random variable, and the histogram of 
all x values is used in order to define an appropriate en- 
semble. For the sake of later discussion we define, be- 
sides the algebraic average ((x)) a also the harmonic av- 
erage as ((x))h = [(1/x)]^ 1 and the geometric average as 
((x)) g = exp[(lna;)]. The result of the resistor network cal- 
culation is labeled as ((x)} (without subscript). 

Our interest is in the circumstances where the strong 
"quantum chaos" assumption of Wigner fails. This would 
be the case if the distribution of x is wide in log scale. 
If x has (say) a log-normal distribution, then it means 
that the typical value of x is much smaller compared with 
the algebraic average. This means that the perturbation 
matrix V nm is effectively sparse {a lot of vanishingly small 
elements). We can characterize the sparsity by the param- 
eter q = ((x}) g /((x)) a . We are going to explain that for 
typical experimental conditions we might encounter sparse 
matrices for which q <C 1 . Then the energy spreading pro- 
cess is similar to a percolation in energy space, and the 
SLRT formula Eq.(3) replaces the Kubo formula. 

Outline. — In what follow we present our model sys- 
tem, analyze it within the framework of SLRT, and then 
introduce an RMT model with log-normal distributed ele- 
ments, that captures the essential ingredients of the prob- 
lem. We show that a generalized resistor network anal- 
ysis for the transitions in energy space leads to a gener- 
alized Variable Range Hopping (VRH) picture (the stan- 
dard VRH picture has been introduced by Mott in [14] and 
later refined by [15] using the resistor network perspective 



of [13]). Our RMT based analytical estimates are veri- 
fied against numerical calculation. Finally wc discuss the 
experimental aspect, and in particular define the physical 
circumstances in which SLRT rather than LRT applies. 
These two theories give results that can differ by orders of 
magnitude. 

Modeling. — Consider a strictly rectangular billiard 
whose eigenstates are labeled by n = (n x ,n y ). The per- 
turbation due to the movement of the 'vertical' wall does 
not couple states that have different mode index n y . Due 
to this selection rule the perturbation matrix is sparse. If 
we deform slightly the potential (Fig. la), or introduce a 
bump (Fig. lb), then states with different mode index are 
mixed. Consequently the numerous zero elements become 
finite but still very tiny in magnitude, which means a very 
wide size distribution featuring a small fraction of large el- 
ements. Similar considerations apply for the circular cav- 
ity of Fig.lc, where an off-center scatterer couples radial 
and angular motion, and which is more suitable for a real 
experiment (but less convenient for numerical analysis). 

Typically the perturbation matrix is not only sparse 
but also textured. This means (see Fig. 2) that there are 
stripes where the matrix elements are larger, and bottle- 
necks where they arc all small. The emergence of tex- 
ture (i.e. non-random arrangement of the sparse large ele- 
ments along the diagonals) is most obvious if we consider 
the geometry of Fig. Id, where wc have a divided cavity 
with a small weakly connected chamber where the driv- 
ing is applied. If the chamber were disconnected, then 
only chamber states with energies E r would be coupled 
by the driving. But due to the connecting corridor there 
is mixing of bulk states with chamber states within en- 
ergy stripes around E r . The coupling between two cavity 
states E n and E m is very small outside of the E r stripes. 
Consequently the near diagonal elements of V nm have wide 
variation, and hence a wide log(x) distribution. 

Coming back to the geometries of Fig.labc, it is some- 
what important in the analysis to distinguish between 
smooth deformation that couples only nearby modes, and 
diffractive deformation that mix all the modes simultane- 
ously: Recalling that different modes have different DOS, 
and that low-DOS modes are sparse within the high-DOS 
modes, we expect a more prominent manifestation of the 
texture in the case of a smooth deformation of a cavity 
that has a large aspect ratio. We later confirm this expec- 
tation in the numerical analysis. 

The SLRT calculation. — As in the standard deriva- 
tion of the Kubo formula, also within the framework of 
SLRT [11,12], the leading mechanism for absorption is as- 
sumed to be FGR transitions. The FGR transition rate is 
proportional to the squared matrix elements |Kim| 2 , and 
to the power spectrum at the frequency to = E n —E m . It 
is convenient to define the normalized spectral function 
F(ui), such that 



RMS(fl) 2 x F(w) 



(4) 
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Contrary to the naive expectation the theory does not lead 
to the Kubo formula. This is because the rate of absorp- 
tion depends crucially on the possibility to make connected 
sequences of transitions. It is implied that both the tex- 
ture and the sparsity of the |V^ m | 2 matrix play a major 
role in the calculation of G. Consequently SLRT leads 
to Eq.(3), where ((...)) is defined using a resistor network 
calculation. Namely, the energy levels are regarded as the 
nodes of a resistor network, and the FGR transition rates 
as the bonds that connect different nodes. Following [12] 
the inverse resistance of a bond is defined as 



2g E 



\V n . 



(E n -E m y 



F(E m —E n ) 



(5) 



and ((\V nm \ 2 )) is defined as the inverse resistivity of the 
network. It is a simple exercise to verify that if all 
the matrix elements are the same, say V,mi| 2 = c, then 
((|Vmn| 2 )} = c too. But if the matrix is sparse or textured 
then typically 



«|V nm | 2 » h « {{\V nm \ 2 )) « «|Krn| 2 ))e 



(6) 



In the case of sparse matrices this is a mathematically 
strict inequality, and we can use a generalized VRH 
scheme which we describe below in order to get an estimate 
for {{x}). If the element-size distribution of log(a;) is not 
too stretched a reasonable approximation is {{x)} « ((x)) g , 
simply because the geometric mean is the typical (median) 
value for the size of the elements. However, if |Ki m | 2 
has either a very stretched element-size distribution, or 
if it has texture, then our VRH analysis below show that 
the geometric average becomes merely an improved lower 
bound for the actual result. 

Analysis. — We consider a particle of mass M in a 
two dimensional box of length L x and width L y , such that 
< x < L x and < y < L y . See Fig. lb. With the driv- 
ing the length of the box becomes R = L x + f(t). The 
Hamiltonian is 



H = diag{£;„} + u{U nm } + f(t){V nm ] 



(7) 



means that F(uo) = 1 up to this cutoff frequency. We 
have also considered (not presented) an exponential line 
shape F(u>) = cxp(— \lo/lu c \), leading to qualitatively sim- 
ilar results. After diagonalization of {E n } + u{U nrn } the 
Hamiltonian takes the form 



H = diag{£„} + f(t){V nm } 



(9) 



where n (not bold) is a running index that counts the 
energies in ascending order. The DOS remains essentially 
the same as for u = 0, namely, 



— ML X L V 
2tt y 



(10) 



The perturbation matrix \V nm \ 2 is sparse and textured 
(see Fig. 2). First we discuss the sparsity, and the effect of 
the texture will be addressed later on. 

Considering first zero deformation [u = 0) it follows 
from Eq.(8) that the non-zero elements of the perturba- 
tion matrix are \V nm \ 2 f» \Mv 2 /L x \ 2 , where v E = ^2E/M. 
The algebraic average of the near diagonal elements equals 
this value (of the large size elements) multiplied by their 
percentage pa. To evaluate po let us consider an energy 
window dE. The number of near-diagonal elements V nm 
within the stripe \E nai , ny -E mxtmy \ < de is gldEde. It is a 
straightforward exercise to find out that the the number of 
non-zero elements (i.e. with n y =m y ) is the same number 
multiplied by pq = [2nMv E Ly] . Consequently 



((\V nm \ 2 ))a 



1 



2ttMv e L v 



L, 



2irL y L 2 x 



(11) 



Somewhat surprisingly this result turns out to be the same 
(disregarding an order unity numerical prcfactor) as for a 
strongly chaotic cavity (see Eq.(I3) of Ref. [3]), as if there 
is no sparsity issue. This implies that irrespective of the 
deformation u, the LRT Kubo result is identical to the 2D 
version of the wall formula (see Sec. 7 of Ref. [3]): 



37T L x 



(12) 



where n = (n a 



is a composite index that labels the Our interest below is not in G LRT but in G SLRT , which 



energy levels E n of a particle in a rectangular box of size 
L x x L y . The deformation is described by a normalized 
Gaussian potential U(x, y) of width (a x , a y ) positioned at 
the central region of the box. Its matrix elements are 
Unm, and it is multiplied in the Hamiltonian by a pa- 
rameter u which signifies the strength of the deformation. 
Note that the limit a — > is well defined and corresponds 
to an "s-scatterer" . The perturbation matrix due to the 
f(t) displacement of the wall is 



(8) 



can differ by many orders of magnitudes. For sufficiently 
small u the large size matrix elements are not affected, 
and therefore the algebraic average stays the same. But 
in the SLRT calculation we care about the small size ma- 
trix elements, that are zero if u = 0. Due to the first-order 
mixing of the levels, the typical overlap |(m|n)| between 
perturbed and unperturbed states is \uU nm /(E n —E m )\. 
The typical size of a small V nm element is the multiplica- 
tion of this overlap (evaluated for nearby levels) by the size 
of a non-zero V nm element. Consequently the small size 
matrix elements are proportional to u 2 . The geometric 
average simply equals their typical size, leading to 



The power spectrum of / is assumed to be constant within 
the frequency range < uj c and zero otherwise. This 



((\Vnm\ 2 )) 9 



M 
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2itL x 
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(13) 
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Motivated by the discussion below Eq.(6) a crude estimate 
for the SLRT result is G SLRT «gx G LRT , where for small 
deformation 



\ 2 )) 9 



«|14 m | 2 »a 



see Eqs.(ll,13) (14) 



It follows from the above (and see Fig. 4) that for 
small deformations q <C 1, and consequently we expect 

LRT 

<G lrt ■ This should be contrasted with the case 
of strongly deformed box for which all the elements are 
of the same order of magnitude and q becomes of order 
unity. Our next task is to further improve the SLRT esti- 
mate using a proper resistor network calculation [d]. 

RMT modeling. — The \V nm \ 2 matrix looks like a 
random matrix with some distribution for the size of the 
elements (see Fig. 3). It might also possess some non- 
trivial texture which we ignore within the RMT frame- 
work. The RMT perspective allows us to derive a quan- 
titative theory for G using a generalized VRH estimate. 
Let us demonstrate the procedure in the case of an ho- 
mogeneous (neither banded nor textured) random matrix 
with log- normal distributed elements. The mean and the 
variance of In (a;) are trivially related to geometric and 
the algebraic averages. Namely, (ln(x)) =ln((a;)) g and 
Var(x) = —2\n(q). Given a hopping range j-EVn — E n \ oj 
we can look for the typical matrix element x^ for con- 
nected sequences of transitions, which we find by solving 
the equation q^ujY^x^) ~ 1, where F(x) is the probability 
to find a matrix element larger than x. This gives 



«z)) g exp 



(15) 



where a = ln(^ E o> c ). From this equation we deduce the 
following: For <?<1, meaning that the distribution is not 
too wide, Xu; « ({x)) g as anticipated. But as the matrix 
gets more sparse (q -C 1), the result deviates from the geo- 
metric average, the latter becoming merely a lower bound. 

The generalized VRH estimate is based on optimization 
of the integral J x^ F(u>) du>. For the rectangular F(ui) 
which has been assumed below Eq.(8) this optimization is 
trivial and gives « x Uo , leading to 



G, 



q exp 



2v/-ln(7 a 



xG T , 



(16) 



where G LRT is given by Eq.(12) and q is given by Eq.(14). 
We have also tested the standard VRH that assumes an 
exponential F(u)) (not presented). 

Numerical results. — The analytical estimates in 
Eqs. (11,13) are supported by the histograms of Fig. 3. For 
each choice of the parameters (AS, a, u) we calculate the 
algebraic, and the geometric and the SLRT resistor net- 
work averages of {|K mi | 2 }. See Fig. 5 and Fig. 6. We also 
compare the actual results for G SLRT with those that were 
obtained from a log-normal RMT ensembles with the same 
algebraic and geometric averages as that of the physical 



matrix [e] . As further discussed in the next paragraph one 
concludes that the agreement of the physical results with 
the associated VRH estimate Eq.(16) is very good when- 
ever the perturbation matrix is not textured, which is in 
fact the typical case for non-extreme aspect ratios. 

In order to figure out whether the result is fully deter- 
mined by the distribution of the elements or else texture 
is important we repeat the calculation for untextured ver- 
sions of the same matrices. The untextured version of a 
matrix is obtained by performing a random permutation 
of its elements along the diagonals. This procedure affects 
neither the bandprofile nor the {ll^ml 2 } distribution, but 
merely removes the texture. In Fig. 5 we see that the phys- 
ical results cannot be distinguished from the untextured 
results, and hence are in agreement with the RMT and 
with the associated VRH estimate. On the other hand, 
in Fig. 6, which is for large aspect ratio, we see that the 
physical results deviate significantly from the untextured 
result. As the width of the Gaussian potential becomes 
larger (smoother deformation), the texture becomes more 
important. These observation are in complete agreement 
with the expectations that were discussed in the modeling 
section. 

Experiment. — As in [4-6] a collection of A~10 6 
atoms, say S5 Rb atoms (M = 1.4 x 10~ 25 fcg), are laser 
cooled to low temperature of T ~ 10/iA^, such that the 
the typical thermal velocity is v E ~ 0.05m/s. The atoms 
are trapped in an optical billiard whose blue-detuned light 
walls confine the atoms by repulsive optical dipolc po- 
tential. The motion of the atoms is limited to the bil- 
liard plane by a strong perpendicular optical standing 
wave. The thickness of the billiard walls (~ 10/j.m) is 
much smaller than its linear size (L ~ 200/jto). The 2D 
mean level spacing is A = g^ 1 ~ 2.5 x 10~ 34 J, which is 
2AHz. One or more of the billiard walls can be vibrated 
with several kHz frequency by modulating the laser in- 
tensity. The dimensionless spectral bandwidth of this 
driving can be set as say w c /A ~ 1000, with an ampli- 
tude ~ lOfim, such that R ~ 0.015m/s. The temperature 
of the trapped atoms can then be measured as a func- 
tion of time by the time-of-flight method. The LRT es- 
timate G LRT ~ 1.3 x 10~ 51 Js/?ti 2 would lead to heating 
rate E ~ 2 x 10~ 27 J/s which is ~ 0.15mA/s. Consider- 
ing (say) the geometry of Fig.lc, the deformation (it) is 
achieved either by introducing an off center optical "spot" , 
or by deforming slightly the optical walls (such precise con- 
trol on the geometry has been demonstrated in previous 
experiments). Having control over u we can have q ~ 10 -5 
that would imply factor 100 suppression, i.e. an estimated 
heating rate of few \iKj 'sec. Such heating rate can be ac- 
curately measured, yielding high sensitivity to the energy 
diffusion process studied here. 

SLRT vs LRT. — Typically the environment intro- 
duces in the dynamics an incoherent relaxation effect. If 
the relaxation rate is strong compared with the rate of 
the externally driven transitions, then the issue of hav- 
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ing "connected sequences of transitions" becomes irrele- 
vant, and the SLRT slowdown of the absorption is not 
expected. In the latter case LRT rather than SLRT is 
applicable. It follows that for finite relaxation rate there 
is a crossover from LRT to SLRT behavior as a function 
of the intensity of the driving. In cold atom experiments 
the relaxation effect can be controlled, and typically it is 
negligible. Hence SLRT rather than LRT behavior should 
be expected. This implies, as discussed above, a much 
smaller absorption rate. Furthermore, as discussed in the 
introduction, one can verify experimentally the signature 
of SLRT: namely, the effect of adding independent driving 
sources is expected to be non-linear with respect to their 
spectral content. 

Conclusions. — In this work we have introduced a 
theory for the calculation of the heating rate of cold atoms 
in vibrating traps. This theory, that treats the diffusion in 
energy space as a resistor network problem, is required if 
the cavity is not strongly chaotic and if the relaxation ef- 
fect is small. The SLRT result, unlike the LRT (Kubo) re- 
sult is extremely sensitive to the sparsity and the textures 
that characterize the perturbation matrix of the driving 
source. For typical geometries the ratio between them is 
determined by the sparsity parameter q as in Eq. (16), 
and hence is roughly proportional to the deformation (it 2 ) 
of the confining potential. If the cavity has a large as- 
pect ratio, and the deformation of the confining potential 
is smooth, then the emerging textures in the perturbation 
matrix of the driving source become important, and then 
the actual SLRT result becomes even smaller. 

By controlling the density of the trapped atoms, or their 
collisional cross section (e.g. via the Feshbach resonance) 
the atomic collision rate can be tuned by many orders 
of magnitude. Their effect on the dynamics can thus be 
made either negligible (as assumed above) or significant, 
thereby serving as an alternative (but formally similar) 
mechanism for weak breakdown of integrability. It follows 
that heating rate experiments can be used not only to 
probe the deformation (u) of the confining potential, but 
also to probe the interactions between the atoms. 
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Fig. 1: Model systems: The atoms are held by a potential 
that may consist of static walls (solid lines), a vibrating wall 
(shaded lines), and bumps (thick points). The numerics has 
been done for (b) with a Gaussian bump. We work with two 
different aspect ratios. For the aspect ratio AS = 20 we take 
L x = 200 and L y — 10. For the aspect ratio AS = 1 we take 
L x = 40 and L y — 40. The position of the Gaussian bump was 
randomly chosen within the region [0.4, 0.6]Li, x [0.4, 0.6]L y . 
The width of the Gaussian is a x — a y — a. We have assumed 
noisy driving with u c = 7A, where A = l/f?E is the mean level 
spacing, and the units were such that M = 1. 
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Fig. 2: Image of the perturbation matrix \V n n 



due to a wall 



displacement of a rectangular-like cavity that has an aspect ra- 
tio AS = 20. The potential floor is deformed due to the pres- 
ence of a cr=0 scatterer with u = 10~ (see text). The matrix 
is both sparse and textured. Lower inset: untextured matrix - 
the elements along each diagonal are randomly permuted. Up- 
per inset: non-sparse matrix with the same band profile - each 
element is generated independently from a normal distribution. 
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Fig. 3: Histograms of matrix elements for different values of 
u for AS = 1 (upper) and AS = 20 (lower) . Here we assume 
a a — scatterer. The vertical lines for u = 10 -2 , 10~ 3 , 10~ 4 
indicate the {{x}) obtained from the LRT algebraic average 
(3 dotted lines that are barely resolved), from the SLRT resis- 
tor network calculation (solid lines), and from the untextured 
calculation (dashed lines) . The geometric mean approximately 
coincides with the peaks, and underestimates the SLRT value 
for the larger AS where the sparsity is much larger. 
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Fig. 4: The sparsity parameter q is plotted versus the 
strength u of the deformation potential for cavities with as- 
pect ratios AS = 1 and AS = 20. We see that for large aspect 
ratio q has some sensitivity to a. As explained in the text 
lrt/Glrt is correlated with q, but for large aspect ratio it 
is even more sensitive to a due to the emergence of textures 
whose presence is not reflected by the value of q. 




Fig. 5: Left panel: The scaled G = {{x}) in the LRT and in 
the SLRT case as a function of u for AS = 1 and different 
smoothness of the deformation. The stars are for the physical 
matrices, while the circles are for their untextured versions (see 
text). The diamonds are for the LRT case. Right panel: The 
SLRT result {{x}) versus the geometric average ({x)) s . These 
are compared with RMT based results, and with the associated 
analytical estimate of Eq.(16). We see the the agreement is 
very good. 




Fig. 6: The same set of plots as in Fig. 5 but for AS = 20. In 
the right panel we clearly see the departure of the physical 
result from the untextured and RMT results, and hence from 
the analytical estimate of Eq.(16). 
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